# This file creates the variables for analyses for parliamentary systems only
# Set working directory
setwd("/Users/ericguntermann/Documents/Papers/Representation of Party Preferences/Replication")

set.seed(123)

load("cses.Rda")
cses <- cses1_4
cses <- cses[-which(cses$countryyear %in% c("Japan 2013")),] # Remove upper house election
## Only keep parliamentary regimes
cses <- cses[!is.na(cses$regime), ]
cses <- cses[cses$regime==0,]
cses <- cses[order(cses$countryyear),]


## Create a matrix of like dislike scores with weights
cses.ld <- cses[,c(4,11:19, 62:64)]

# Remove people who expressed no preferences
cses.ld[is.na(cses.ld)] <- 11
rowtot <- rowSums(cses.ld[,2:10], na.rm=T)
cses.ld <- cses.ld[!rowtot==99,]
cses.ld[cses.ld==11] <-NA
cses <- cses[!rowtot==99,]

## Create overall weight (sample, demographic, and political)
cses.ld$weight <- cses.ld$sampleweight*cses.ld$demoweight*cses.ld$polweight

# Criterion 1: What proportion of people in each country's first choice party ended up in government?
cses.ld2 <- data.frame(cses.ld)
cses.ld2[is.na(cses.ld2)] <- 0


# Create a vector of first choice parties for all individuals
# I used which.is.max, which breaks ties at random
library(nnet)
first <- apply(cses.ld2[,2:10],1, which.is.max)

# Cabinet data
incabinet <- cses[,c(4,21:29)]


# Is the most liked party in government (for each respondent)?
firstincabinet <- vector()
for(i in 1:length(first)){
  firstincabinet[i] <- incabinet[i,first[i]+1]
}

# Get proportion whose most liked party is in cabinet for each election
propfic <- vector()
for (i in 1:length(unique(cses$countryyear))){
  propfic[i] <- weighted.mean(as.numeric(firstincabinet[cses$countryyear == unique(cses$countryyear)[i]]), cses.ld$weight[cses.ld$countryyear==unique(cses.ld$countryyear)[i]], na.rm=T)
}


# Criterion 2: 
## Calculate mean like-dislike scores by country-year

cses.ld.ag <- matrix(nrow=length(unique(cses.ld$countryyear)), ncol=9)
for(i in 1:length(unique(cses.ld$countryyear))){
  for(j in 1:9){
    cses.ld.ag[i,j] <- weighted.mean(cses.ld[c(cses.ld$countryyear==unique(cses.ld$countryyear)[i]),1+j], cses.ld$weight[cses.ld$countryyear==unique(cses.ld$countryyear)[i]], na.rm=T)
  }
}

cses.ld.ag <- as.data.frame(cses.ld.ag)


## Create a vector of most liked parties (aggregate)

cses.ld.ag[,10] <- apply(cses.ld.ag[,1:9], 1, which.max)

mostliked <- cses.ld.ag$V10


## Create aggregate cabinet data matrix

incabinetag <- data.frame(matrix(nrow=length(unique(cses$countryyear)), ncol=9))
for (i in 1:length(unique(incabinet$countryyear))){
  for (j in 1:9){
    incabinetag[i,j] <- incabinet[c(incabinet$countryyear==unique(incabinet$countryyear)[i]),j+1][1]
  }
}
colnames(incabinetag) <- c("incabinet1", "incabinet2", "incabinet3", "incabinet4", "incabinet5", "incabinet6", "incabinet7", "incabinet8", "incabinet9")



## Create logical vector indicating whether the most preferred party is in cabinet
for (i in 1:length(unique(cses$countryyear))){ 
  incabinetag$isin[i] <- incabinetag[i, mostliked[i]]==1
}


# Criterion 3: How liked are governing parties compared to opposition parties?
## Create matrix of governing parties
gov <- cses[, c(4, 21:29)]

## Create matrix of non-governing parties
nongov <- cses[, c(4, 21:29)]
nongov[nongov==1] <- 2
nongov[nongov==0] <- 1
nongov[nongov==2] <- 0

## Create matrix of aggregate weights for governing parties
gweights <- cses[,50:58]*gov[,2:10]
gweightstot <- apply(gweights, 1, sum, na.rm=T)
gweights <- gweights/gweightstot


## Create matrix of weights for non-governing parties
ngweights <- cses[, 50:58]*nongov[, 2:10]
ngweightstot <- apply(ngweights, 1, sum, na.rm=T)
ngweights <- ngweights/ngweightstot


## Overall like dislike score for parties in parliament that are in cabinet
cses.ld2 <- cses.ld[, 2:10]
ldcabinet <- gweights*cses.ld2
ldoverallc <- apply(ldcabinet, 1, sum, na.rm=T)

## Overall like dislike score for parties in parliament that aren't in cabinet

ldnoncabinet <- ngweights*cses.ld2
ldoverallnc <- apply(ldnoncabinet, 1, sum, na.rm=T)

## Create vector of differences between like-dislike scores for parties in cabinet and parties not in cabinet

ldcnc <- data.frame(ldoverallc, ldoverallnc)
lddiff <- vector()
for (i in 1:length(cses$countryyear)){
  lddiff[i] <- ldcnc[i,1]-ldcnc[i, 2]
}

lddiffag <- vector()
for(i in 1:length(unique(cses$countryyear))){
  lddiffag[i] <- weighted.mean(lddiff[cses$countryyear==unique(cses$countryyear)[i]], cses.ld$weight[cses.ld$countryyear==unique(cses.ld$countryyear)[i]], na.rm=T) 
}
lddiff <- lddiffag



# Calculate number of parties in government

ningov <- function(x){
  return(sum(x==1,na.rm=T))
}
ngov <- apply(incabinetag[,1:9], 1, ningov)


## Show data
# Create dataset
criteria <- data.frame(countryyear=unique(cses$countryyear), mostliked=incabinetag$isin, lddiff, propfic)
for (i in 1:length(unique(cses$countryyear))){
  criteria$proportional[i] <- cses$proportional[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$regime[i] <-  cses$regime[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$mdm[i] <-  cses$mdm[cses$countryyear==unique(cses$countryyear)[i]][1]
}
for (i in 1:length(unique(cses$countryyear))){
  criteria$gallagher[i] <- cses$gallagher[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$freedomhouse[i] <-  cses$freedomhouse[cses$countryyear==unique(cses$countryyear)[i]][1]
}

for (i in 1:length(unique(cses$countryyear))){
  criteria$gdppercap[i] <-  cses$gdppercap[cses$countryyear==unique(cses$countryyear)[i]][1]
}


criteria$ngov <- ngov
criteria$coalition <- as.factor(ngov>1)
# Adjust number of parties in Australia. The Liberal and National Parties ran together.
criteria$ngov[criteria$countryyear=="Australia 1996"] <- 1
criteria$ngov[criteria$countryyear=="Australia 2004"] <- 1
library(foreign)
write.dta(criteria, "criteria_parl.dta")


library(xtable)



# This code produces table 11 of online appendix

dat <- read.dta("criteria_parl.dta")
criteria <- dat
criteria2 <- dat[-which(dat$countryyear %in% c("Thailand 2001", "Thailand 2011","Japan 1996", "Japan 2013")),]


#Criterion 1

tab1 <- data.frame(round(mean(criteria2$propfic[criteria2$proportional==0])*100,1), round(mean(criteria2$propfic[criteria2$proportional==1])*100,1), round(mean(criteria$propfic)*100,1))
colnames(tab1) <- c("Non-PR", "PR", "Overall")
rownames(tab1) <- c("Proportion most liked in cabinet (%)")

#Criterion 2
tab2 <- data.frame(t(round(prop.table(with(criteria2, table(mostliked==FALSE, proportional)),2)*100,1)[1,]), round(prop.table(table(criteria$mostliked==TRUE))[2]*100,1)) 
colnames(tab2) <-c("Non-PR", "PR", "Overall")

#Criterion 3
tab3 <- data.frame(round(mean(criteria2$lddiff[criteria2$proportional==0]),2), round(mean(criteria2$lddiff[criteria2$proportional==1]),2), round(mean(criteria$lddiff),2))
colnames(tab3) <- c("Non-PR", "PR", "Overall")
rownames(tab3) <- c("Evaluation of governing vs. opposition parties")


taball <- as.matrix(rbind(tab1, tab2, tab3))
rownames(taball)[2] <- "Most liked party overall in cabinet"

xtab <- xtable(taball, caption="Performance of PR and non-PR Systems on Three Criteria")
sink("taballparl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), only.contents=T)
sink(NULL)

# Tables 12 and 13


## Criterion 1 (proportional)
criteria <- read.dta("criteria_parl.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d1 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))

m1 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p1 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse")

i1 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}


## Criterion 1 (gallagher)
library(foreign)
dat <- read.dta("criteria_parl.dta")
d2 <- with(dat, list(propfic=propfic, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(dat$countryyear))
m2 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i2 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p2 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 1 (mdm)
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d3 <- with(criteria, list(propfic=propfic, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d3$N <- length(unique(criteria$countryyear))

m3 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p3 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i3 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}

## Criterion 1 (proportional) with number of parties
criteria <- read.dta("criteria_parl.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d4 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))

m4 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]+ b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p4 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse", "b.ngov")

i4 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}


d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=20000, n.thin=10, n.chains=2)

results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])


# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion1parl.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT
#Model 1 (Proportional)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co1 <- c(mean(b), mean(b.proportional), NA, NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), NA, NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), NA, NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co2 <- c(mean(b), NA, mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 3 (MDM)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co3 <- c(mean(b), NA, NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P","Mean", "SD","P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1parl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

#Model 4 (Proportional with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]

co1 <- c(mean(b), mean(b.proportional), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co1,sd1,p1)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Number of parties", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1bparl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)


# Tables 14 and 15

## Criterion 2 (gallagher)
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d1 <- with(dat, list(mostliked=mostliked, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(dat$countryyear))
m1 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i1 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p1 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 2 (mdm)
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d2 <- with(criteria, list(mostliked=mostliked, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(criteria$countryyear))

m2 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p2 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i2 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}



## Criterion 2 (gallagher) with number of parties
library(foreign)
dat <- read.dta("criteria_parl.dta")
d3 <- with(dat, list(mostliked=mostliked, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d3$N <- length(unique(dat$countryyear))
m3 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i3 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}
p3 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse", "b.ngov")



## Criterion 2 (mdm) with number of parties
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d4 <- with(criteria, list(mostliked=mostliked, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))
m4 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p4 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse", "b.ngov")


i4 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}

d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=20000, n.thin=10, n.chains=2)


results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])

# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion2parl.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT


#Model 1 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co1 <- c(mean(b), mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b),sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (MDM)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co2 <- c(mean(b), NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))



library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2)                                                                                                                
rownames(model.effects) <- c("Intercept", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 2", digits=c(2))
sink("criterion2parl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

#Model 3 (Gallagher with ngov)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co2 <- c(mean(b), mean(b.gallagher), NA, mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), sd(b.gallagher), NA, sd(b.ngov),  sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 4 (MDM with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co3 <- c(mean(b), NA, mean(b.mdm), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, sd(b.mdm), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)),ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Gallagher*", "MDM*", "Parties in Govt", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 2", digits=c(2))
sink("criterion2bparl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

# Tables 16 and 17

## Criterion 3 (proportional)
criteria <- read.dta("criteria_parl.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011","Japan 1996", "Japan 2013")),]
d1 <- with(criteria, list(lddiff=lddiff, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))

m1 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p1 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse")

i1 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}



## Criterion 3 (gallagher)
library(foreign)
dat <- read.dta("criteria_parl.dta")
d2 <- with(dat, list(lddiff=lddiff, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(dat$countryyear))
m2 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i2 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p2 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 3 (mdm)
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d3 <- with(criteria, list(lddiff=lddiff, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d3$N <- length(unique(criteria$countryyear))

m3 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p3 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i3 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}

## Criterion 3 (mdm) with number of parties
library(foreign)
criteria <- read.dta("criteria_parl.dta")
d4 <- with(criteria, list(lddiff=lddiff, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))
m4 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p4 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse", "b.ngov")


i4 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}


d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=20000, n.thin=10, n.chains=2)


results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])

# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion3parl.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT
#Model 1 (Proportional)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co1 <- c(mean(b), mean(b.proportional), NA, NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), NA, NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), NA, NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co2 <- c(mean(b), NA, mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 3 (MDM)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co3 <- c(mean(b), NA, NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P","Mean", "SD","P")
xtab <- xtable(model.effects, caption="Criterion 3", digits=c(2))
sink("criterion3parl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)


#Model 4 (MDM with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co1 <- c(mean(b), mean(b.mdm), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.mdm), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)),ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))



library(xtable)
model.effects <- data.frame(co1,sd1,p1)                                                                                                                
rownames(model.effects) <- c("Intercept", "MDM*", "Number of parties", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 3", digits=c(2))
sink("criterion3bparl.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)




